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The conductivity of highly charged membranes is nearly constant, due to counter-ions screening 
pore surfaces. Weakly charged porous media, or "leaky membranes", also contain a significant 
concentration of co-ions, whose depletion at high current leads to ion concentration polarization 
and conductivity shock waves. To describe these nonlinear phenomena the absence of electro- 
osmotic flow, a simple Leaky Membrane Model is formulated, based on macroscopic electroneutrality 
and Nernst-Planck ionic fluxes. The model is solved in cases of unsupported binary electrolytes: 
steady conduction from a reservoir to a cation-selective surface, transient response to a current 
step, steady conduction to a flow-through porous electrode, and steady conduction between cation- 
selective surfaces in cross flow. The last problem is motivated by separations in leaky membranes, 
such as shock electrodialysis. The article begins with a tribute to Neal Amundson, whose pioneering 
work on shock waves in chromatography involved similar mathematics. 
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Dedication by Martin Z. Bazant 

This article is dedicated to the memory of Neal R. 
Amundson, the "father of modern chemical engineer- 
ing" pQ, who brought mathematical rigor to the fields 
of transport phenomena and reactor engineering [51 [3J. 
His education, teaching and research were truly interdis- 
ciplinary, long before that term came into fashion. His 
early education (BS 1937, MS 1941) was in Chemical En- 
gineering, the field of his primary faculty appointments 
at University of Minnesota (1949-1977) and University of 
Houston (1977-2011) and lifelong professional focus. He 
is famous for leading the Department of Chemical Engi- 
neering at Minnesota to lasting national prominence, as 
its chair for 25 years - starting at age 33, only two years 
after being hired. It is perhaps surprising then, that his 
most advanced degree at Minnesota (PhD 1945) and his 
early teaching as an Assistant Professor (1945-1947) were 
not in Chemical Engineering, but in Mathematics. 

Amundson revolutionized the way that chemical en- 
gineers design systems for separations, heat transfer 
and adsorption, by replacing empirical principles with 
rigorous models based on partial differential equations 
(PDE). His PhD thesis (1945) [4] and early papers 
(1948-1952) [5HZ] involved analytical solutions of PDEs 
for flow-through adsorption in porous media, which are 
still used today in chromatography, electrophoresis and 
ion exchange. His work built upon the seminal paper 
of Thomas (1944) [5] and preceded those of Goldstein 
(1953) [Hl[in]i which are better known in applied math- 
ematics [TT] , 

A major achievement of Amundson was to show that 
flow-through adsorption processes described by the PDE, 
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where c is the flowing concentration and /(c) is total 



(adsorbed + flowing) concentration per volume in local 
equilibrium, lead to nonlinear kinematic waves [3J lllj . 
Neglecting diffusion (D = 0), the model reduces to a first- 
order quasi-linear PDE, which can be solved in implicit 
form, 
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by Lagrange's Method of Characteristics, as explained in 
a series of papers by Rhee, Aris and Amundson (1970- 
1972) [L2HT4] . For most boundary and initial condi- 
tions, kinematic waves eventually "break" and produce 
shock waves, or propagating discontinuities, which are 
sharp in the limit D = 0. Lapidus and Amundson 
(1952) 7 showed that these discontinuities are broad- 
ened by diffusion (D > 0) around the same time that 
Hopf [TS] and Cole [TB] famously solved Burgers equation 
in fluid mechanics [11] . Amundson, Aris and Swanson 
(1965) [17] rigorously explained and analyzed the sharp 
concentration bands arising in chromatography and ion 
exchange as concentration shocks in porous exchange 

beds Ha using. 

Upon being invited to contribute to this special issue, 
I set out to learn more about the man and his life's work. 
At first, I was struck by the unusual parallels with my 
own career, since I too began graduate study and teach- 
ing in mathematics before moving to chemical engineer- 
ing and holding a joint appointment. What impressed me 
most, however, were the parallels in research. Without 
knowing Admunson's seminal work on exchange beds, A. 
Mani and I recently developed a theory of "deionization 
shocks" in porous media [19] that bears intriguing math- 
ematical similarities; surface conduction within the pores 
and bulk ionic current play the roles of surface adsorption 
and pressure-driven flow, respectively. This phenomenon 
was first discovered in microfluidics by Mani, Zangle and 



Santiago (2009) [20H23"] . but its extension to porous me- 
dia - analogous to Amundson's work - is more general 
(e.g. decoupling the directions of flow and current) and 
paves the way for practical applications. Indeed, the first 
experiment in my new laboratory in the Department of 
Chemical Engineering [24] applies the shock theory to 
water deionization by "shock electrodialysis" [25] and re- 
lies on mathematical analysis for finite-length pores [2(5] 
to interpret the data, as elaborated in this article. 

Since Amundson and I both began our careers in math- 
ematics, it is perhaps not surprising that we ended up 
doing similar research in chemical engineering. At MIT, 
I used to teach 18.311 Principles of Applied Mathemat- 
ics, which introduced PDEs to undergraduates, starting 
with first-order quasilinear equations (applied to traffic 
flow) |27j . In contrast, chemical engineering courses up 
to the graduate level (such as 10.50 Analysis of Trans- 
port Phenomena, which I teach with W. Deen 28J follow- 
ing Amundson's co-teaching model [Tj) mainly cover the 
solution of linear parabolic PDEs, including the Finite 
Fourier Transform method championed by Amundson [2] . 
As a result, most chemical engineering students today do 
not know the Method of Characteristics, even though it 
is the mathematical basis for theories of chromatography 
[3 D3 > gas dynamics [TTJ [25] , electrokinetic soil remedi- 
ation 30-33] . capillary electrophoresis [341 - I38] . and ion 
concentration polarization in microchanncls [20 22] and 
porous media [HI I3S] (the focus of this article) . 

Perhaps this special issue of AIChE Journal can serve 
as a call to reinvigorate the teaching of mathematical 
methods introduced by Amundson to our field, which 
provide physical insights and useful formulae, too often 
overlooked in the Computer Age. □ 



Introduction 

When current is passed through an electrolyte to an 
ion selective surface (such as an ion-exchange membrane, 
micro/nanochannel junction or electrode) , the passage of 
certain ions, and the rejection of others, generally lead to 
concentration variations and voltage losses (or internal 
resistance), known as "ion concentration polarization" 
(ICP). Under classical assumptions of electroneutrality 
without convection or homogeneous reactions, the cur- 
rent is limited by electrodiffusion, when the concentration 
of the active species vanishes at the selective surface [40] . 
In a neutral binary electrolyte, the current appears to be 
limited by diffusion alone, since the concentration pro- 
files evolve according to a pure diffusion equation, but 
electromigration and diffusion conspire to determine the 
effective (ambipolar) diffusivity [5T] . Both species diffuse 
in the same direction, but electromigration enhances the 
flux of the active species and opposes the flux of the in- 
active species (and cancels it in steady state). 

Despite this theoretical speed limit, over limiting cur- 
rent (OLC) has been observed in a variety of sys- 



tems involving membranes, porous media, and mi- 
cro/nanochannels. Elucidating mechanisms for OLC re- 
mains a central question in membrane science and chem- 
ical engineering [12]. In free electrolyte solutions, there 
are two fundamental mechanisms for OLC - chemical and 
physical - each of which affects ICP in different ways. 

Chemical mechanisms for OLC involve the production 
of additional ions (from solvent decomposition, H + and 
OH~ in water) and/or the loss of surface selectivity (from 
charge regulation of a membrane or nanochannel or from 
side reactions at an electrode) in order to reduce ICP and 
maintain ionic conductivity at high currents [42H44] . An- 
dersen et al. (2012) [45] recently showed that both phe- 
nomena are needed to achieve significant OLC via the 
phenomenon of "current-induced membrane discharge" 
(CIMD). In particular, for aqueous systems, CIMD can 
result from bulk water splitting coupled to charge reg- 
ulation of the membrane, e.g. by proton adsorption in 
anion exchange membranes. 

Physical mechanisms for OLC involve the amplifi- 
cation of a different transport mechanism, which al- 
lows ions to reach the selective surface faster than by 
quasi-neutral electrodiffusion in the region of strongest 
ICP near the selective surface. The best-known exam- 
ple is the Rubinstcin-Zaltzman electro-osmotic instabil- 
ity 46-49J, which has recently been observed in mi- 
cro/nanofluidic systems [50H52] . In porous media and 
microchannels, the presence of charged double layers on 
the side walls, aligned with the direction of current, al- 
lows OLC to be sustained by the additional transport 
mechanisms of surface conduction (SC) [26] and electro- 
osmotic flow (EOF) [251 E3 ■ At lower (under-limiting) 
currents, ICP has also been observed experimentally in 
an electroosmotic pump with a porous glass frit [54] and 
in porous electrodes in a system designed for capacitive 
desalination 55 . 

Advances in microfluidics over the past ten years 
have enabled the direct observation of ICP during 
OLC. Steady, sharply defined depletion zones, some- 
times containing internal electro-osmotic vortices, have 
been observed near the junctions of microchannels and 
nanochannels by by J. Han's group since 2005 [56-58 
and have been shown to be affected by the microchannel 
geometry [5U[SD]. Of particular note for this work, Mani, 
Zangle and Santiago (2009) have shown that in very thin 
(1 //m) channels, these depletion interfaces can propagate 
as shock waves under constant current [191121] . Zangle et 
al. (2010) 22] give an insightful review of these experi- 
ments and the theory behind them, which is mathemat- 
ically similar to Amundson's theory of chromatography, 
as noted above. 

While the original experiments and theory were limited 
to single microchannels, Mani and Bazant (2011) 19J pre- 
dicted the possibility of propagating deionization shocks 
in porous media and formulated general nonlinear PDEs 
to describe volume-averaged ICP at the macroscopic 
scale, driven by SC within charged nanopores, neglect- 
ing EOF as a first approximation. They obtained an- 



alytical similarity solutions for power-law variations in 
microstructure and analyzed the internal structure and 
dynamical stability of deionization shocks in three di- 
mensions. The formal volume averaging of their macro- 
scopic ICP model is analogous to Amundson's theory of 
surface adsorption in fixed beds [5HZ1 H21 HZ| and Hclf- 
ferich's early models of ion exchange membranes |61) . but 
the PDEs are solved under general conditions of strong 
ICP with diffusion. Although porous ICP equations also 
provide simple area-averaged descriptions of microfluidic 
devices, they are also more general because the flow and 
current directions can be decoupled in three dimensions, 
opening some new possibilities for separations. 

In this article, we develop a fundamental picture of 
ICP dynamics in finite-size porous domains, including 
effects of simultaneous pressure-driven fluid flow and ap- 
plied current, leading to two-dimensional concentration 
variations. The analysis is based on the PDEs of Mani 
and Bazant |19j for nanopores (or nanochannels) in the 
SC regime with forced convection, neglecting electro- 
osmotic flows that dominate in larger pores at the mi- 
cron scale 26J. Borrowing a term of A. Yaroshchuk [B"2"] . 
we will refer to this as the "Leaky Membrane Model". 
All of the example calculations here are motivated by 
experiments in our group aimed at establishing surface- 
transport mechanisms for OLC and harnessing ICP dy- 
namics in porous media for novel separations [24 . 

Before we begin, we would like to draw attention to 
the recent work of Yaroshchuk, connecting these ideas to 
classical membrane science [ST] via in a theory of "leaky 
membranes" and performing some similar transient [39] 
and steady state [621 163) calculations in one dimension, 
without flow. His analysis is based on "virtual concentra- 
tions" in local thermodynamic equilibrium with a hypo- 
thetical ionic reservoir across each thin slice of the porous 
medium. The results can be left in general form or con- 
nected to specific quasi-cquilibrium local models, such as 
the Poisson-Boltzmann model with fixed surface charge 
in a straight channel, with thin or thick double layers [B"4"] . 
This is an analytical limit of the full model of Mani et 
al 20J for thick double layers with effective longitudinal 
transport coefficients obtained numerically by integrating 
over the cross section. Our approach uses the physical 
volume-averaged concentration variables (defined in the 
next section) and the slowly varying, macroscopic part 
of the potential of mean force. As such, interfacial volt- 
ages at the ends of the porous domain must be added 
to describe experimental data, but this can be done ac- 
curately by modeling or measuring the electrochemical 
series resistances and open circuit voltage [2~4l . 



Leaky Membrane Model 

General Formulation 



charge/area a s filled with an electrolyte, containing ions 
of charge z-^e and concentration Cj (per pore volume). An 
example of a porous silica glass frit is shown in Fig. [T] 
Charge conservation implies, 
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which can be written as a balance of surface charge per 
pore volume, 
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with the electrolyte charge density. 

The electroneutrality condition Q determines the rel- 
ative importance of surface charge. Let cq be a typical 
concentration of co-ions (of the same sign as the surface 
charge), which sets the scale for neutral salt permeating 
the porous medium. There are two limiting cases, illus- 
trated in Fig. flTc) . In the membrane limit, \p s \ ^$> eco, 
co-ions are strongly excluded, and the porous medium 
maintains a large, constant conductivity from nearly uni- 
formly distributed counter-ions (of the opposite sign as 
the surface charge). In the opposite limit, |p s | <C eco, 
naively one would expect classical electrodiffusion of a 
quasi-neutral electrolyte, only with diffusivities rescaled 
to account for the tortuosity of the medium. This is 
indeed the case at low currents, but it turns out that 
the surface charge is a singular perturbation. As we 
shall see, a "leaky membrane" with \p s \ > 0, no matter 
how small, is fundamentally different from an uncharged 
porous medium with p s = 0. 

The simplest approximation for the nonlinear dynam- 
ics of the electrolyte in a charged porous medium, which 
we call the Leaky Membrane Model (LMM), combines 
the macroscopic electroneutrality condition Q with ho- 
mogenized Nernst-Planck equations [19] 
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where F^ is the macroscopic flux density of species i (per 
cross-sectional pore area), Di is the macroscopic chemi- 
cal diffusivity (with tortuosity correction |65j). (f> is the 
slowly varying, non-equilibrium part of the electrostatic 
potential of mean force, u is the mean fluid velocity, 
and R t is the mean reaction rate (per volume) producing 
species i. The current density is 



^2 z i eF i 



(7) 



Consider a charged porous medium of porosity e p , 
internal pore surface area/ volume a p , pore surface 



(per cross-sectional pore area). For binary electrolytes, 
the LMM can also be cast in terms of bulk and surface 
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small pore, thick double layers 



large pore, thin double layers 
bulk conductivity K b ~ c_ 
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FIG. 1: Physical picture of a leaky membrane, (a) SEM image of a silica glass frit with 557 nm mean pore size, which can 
sustain over-limiting current and deionization shocks (Deng et al. |24).'l (b) Fixed surface charges (red) of density a s per pore 
area and p s per total volume, and mobile counter-ions (orange) and co-ions (yellow) in the pores of concentration c± per total 
volume, (c) Sketch of the quasi-equilibrium ion profiles in small and large pores, relative to the Debye length Ad. The surface 
conductivity scales with the total screening charge (orange areas), and the bulk, depletable conductivity scales with the co-ion 
charge (yellow areas). 



conductivities (defined below) [T9"] . 

In a concentrated solution, the chemical diffusivity de- 
pends on the ionic concentrations [HJ |S3] , 



D, = D ( 1 
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(8) 



where D® is the tracer diffusivity in a dilute solution 
(which also generally depends on concentration [6"5T - 
[67] ) and ji is the molar activity coefficient. Using ^ 
in (p)]) leads to "modified Poisson-Nernst-Planck equa- 
tions" [5H1 [SH], which are used to account for ther- 
modynamic effects in the nonlinear dynamics of elec- 
trolytes [70l [71] . In a general formulation based on non- 
equilibrium thermodynamics [66j . the ionic fluxes and 
reaction rate are related to electrochemical potentials, 
defined as variational derivatives of the total free energy 
functional. For systems with phase transformations (such 
as precipitation or phase separation), the Nernst-Planck 
equations ^ become generalized Cahn-Hilliard-reaction 
models [66] coupled with macroscopic quasi-neutrality 
Q in the LMM framework. Coupled diffusive fluxes, e.g. 
friction between different species in a Stefan-Maxwell for- 
mulation, can also be important for large ions or charge 
colloids in porous media [72] . 

Here we focus mostly on dilute solutions, where Di is 
constant or simply a function of the local salt concentra- 
tion. 



Theoretical Justification 

In the absence of flow and reactions, the LMM can 
be derived from the microscopic Poisson-Nernst-Planck 
(PNP) equations within the pores by taking the limit 
of thin double layers and area averaging [20] , by formal 



homogenization of the microscopic PDEs for arbitrary 
double layer thickness [73], and by assuming local ther- 
modynamic quasi-equilibrium at the pore scale [39] . The 
microscopic potential of mean force within the pores is 
4> + V'eq, where <j> is the slowly varying part reflecting 
macroscopic departures from equilibrium and ip eq is the 
rapidly varying correction due to quasi-equilibrium lo- 
cal interactions between the ions and the surface charge, 
constrained by the slowly varying mean ion concentra- 
tions Ci. In the microscopic PNP equations, the surface 
charge per internal pore area, a s , enters via the electro- 
static boundary condition on the pore walls, but after ho- 
mogenization the macroscopic Nernst-Planck equations 
( 5 ) are simply augmented by a quasi-neutrality condition 
(4) that includes the surface charge per volume p s . Ef- 
fectively, the local quasi equilibrium charge fluctuations 
associated with %\) eq are "integrated out" by homogeniza- 
tion to macroscopic length scales in the porous medium. 

Fluid flow in charged porous media is much more com- 
plicated to homogenize rigorously, and no simple ap- 
proximations are currently available. The difficulty is 
that flows within the pores are strongly coupled to the 
ion profiles via (locally linear) electrokinetic phenomena 
and lead to complex dispersion effects, modifying Di by 
nonuniform convection in the porous medium. Classi- 
cal Taylor dispersion [53] is often negligible compared to 
internal electro-osmotic convection [26] and eddy disper- 
sion [24] . In the simplest version of the LMM considered 
here, we neglect electroconvection and dispersion in ^ 
and simply assume an imposed pressure-driven flow. 

In this work, we also neglect the reaction rate Ri, 
but reactions are important in many situations, such as 
electrokinetic remediation [30H33] and porous electrode 
charging [65] [111 [75], and provide an additional source 
of nonlinearity. In particular, the surface charge density, 
<7 S , is generally a function of the local electrolyte com- 



position via the specific adsorption of ions. This phe- 
nomenon of "charge regulation" is crucial for the quanti- 
tative interpretation of shock electrodialysis experiments 
using LMM [21] and also underlies the theory of current- 
induced membrane discharge 45 . Assuming fast adsorp- 
tion kinetics, the LMM then closely resembles Amund- 
son's classical models of exchange beds of the form (IT]), 
except that the LMM also accounts for the electroki- 
netic coupling between ion adsorption via the macro- 
scopic charge balance Q. 



reference solution of salt concentration c$. As discussed 
above, the key parameter is p s , which determines to what 
extent the porous medium acts like a "good membrane" 
with high conductivity and selectivity for counter-ions 
(\p s \ 3> 1)- We are mainly interested in "leaky mem- 
branes" with < \p s \ -C 1, which become depleted at 
high currents, leading to complex nonlinear dynamics. 



Surface Conduction in a Leaky Membrane 



Uniform Membrane Charge in a Binary Electrolyte 

In order to highlight the nonlinear dynamics of ion 
transport in porous media, we adopt the simplest form of 
the LMM. As noted above, we neglect reactions (Ri = 0) 
and charge regulation (-£*■ = 0). We also impose a 
pressure-driven Darcy flow u without accounting for dis- 
persion or electrokinetic phenomena. In particular, we 
consider only the representative cases of uniform flow, 
either parallel or perpendicular to the applied current. 
We further assume a uniform porous medium with con- 
stant microstructure and charge {p s ,t p =constants). For 
ease of calculations, we also make the standard theoret- 
ical assumption of a symmetric binary z : z electrolyte 
with equal ionic diffusivities D (including the tortuosity 
correction). See Mani and Bazant [TS] for extensions to 
asymmetric binary electrolytes and nonuniform porous 
media (with uniform flow) and Andersen et al. |45j for 
the full LMM (without flow) for a multicomponent (four 
species) electrolyte in a leaky membrane with charge reg- 
ulation (proton adsorption) and homogeneous reactions 
(water self- ionization) . 

With these assumptions, the LMM takes the form 
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in terms of the dimensionless variables, x = x/L, V = 
LV, I = tD/L 2 , c± = c ± /c , 4> = ze(f>/k B T, F ± = 
F±L/Dco and u = u/U, for a geometrical length scale L 
and characteristic velocity U. There are two dimension- 
less groups that control the solution, the Peclet number 
(ratio of convection to diffusion) 



Pc = 



UL 
~D 



and the dimensionless charge density, 

a r ,(j. 
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(13) 



which is the ratio of fixed surface charges to mobile ionic 
charges per volume, if the pores were filled with a neutral 



Before we proceed to the analysis, we comment on 
the non-standard definition of "surface conduction" in 
the LMM 19J. In classical electrokinetic systems, the 
term 'surface conduction' refers to the excess conduction 
(from electromigration and electro-osmotic convection) 
that arises from increased ion concentrations in the elec- 
tric double layer (EDL) [751 EZ] • Surface conduction in 
this case can be found by taking the total conduction 
and subtracting the conduction that would be found in 
the absence of an EDL. This definition has a long history 
in electrokinetics from pioneering work of Smoluchowski, 
Bikcrman and Dukhin|78 84 . For a particle or pore of 
characteristic length a, the Dukhin number 
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CLKb 



(14) 



is the ratio of excess surface conductivity, k s , to bulk 
conductivity, Kb- 

While surface conduction is traditionally defined in 
terms of the excess ion concentration EDL, the more rel- 
evant definition for a leaky membrane is in terms of the 
total surface charge density, as shown in Fig. [T] The dif- 
ference between these two values is shown in Figure [2] 
Let r± be the total excess surface concentration (Fig- 
ured). This can be written as 



T± = / ( c ± - c bu i k )dx 

for a binary, univalent electrolyte, where c_, c+, and 
Cbuik are the negative, positive, and bulk ion concentra- 
tions, respectively, and x is the distance from the charged 
surface. For a negatively charged surface, T + > and 
r_ < 0. In classical electrokinetics, the excess surface 
conductivity will come from the excess neutral salt con- 
centration, w [BID |5B] • 

w = (c+ + C- - 2c bulk )dx = r + + r_ 

In a leaky membrane, however, the total surface ion con- 
centration, or double-layer charge density, q, plays an 
important role and is given by 

q = / (c+ - c-)dx = T + - T_ = -cr s 

Here we consider an EDL at equilibrium and examine 
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FIG. 2: Ion Distribution Near a Charged Surface, a) Excess ion concentration, b) Excess neutral concentration, c) Surface 
charge density 



a different mechanism for conduction along a charged 
surface. This surface conduction is in addition to, rather 
than in excess of, the conduction through the neutral 
bulk electrolyte. Throughout this paper, only the surface 
conduction due to the total surface charge density will be 
discussed, and will be referred to as SC. 

The importance of surface phenomena generally in- 
creases with the surface to volume ratio. In the case 
of SC, the role of surface charge is controlled by p s , the 
ratio of surface charge to bulk ionic charge per volume, 
which becomes non- negligible in submicron pores, espe- 
cially at low electrolyte concentrations. It is important 
to note that this dimensionless group, while similar, is 
not the same as the Dukhin number 19 . The Dukhin 
number depends on w, while p~ s depends on d s = — q. 
It is possible for w to be very small while maintaining a 
large q value, resulting in a large value for p s and a small 
Dukhin number. 



Uniform Current without Flow 



Eqs. (l9|-(ll), into a dimensionless PDE 

d 2 d> 



dc 
dt 
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(15) 



for the depletable salt concentration c = £- (which, as 
explained above, is equal to the co-ion concentration) and 
a constraint for the uniform, time- varying current, 



I = - (C - Ps) 



d.f 



(16) 



obtained by integrating the cation transport equation. 
The dimensionless current 



/ = 



IL 
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(17) 



is carried by cations and scaled to the limiting current 
for the case of an "ideally leaky" membrane, p s = 0. 
The boundary conditions fix the reservoir concentration, 
c(0, t) = 1, and impose zero anion electro-diffusive flux 
at the cation-selective surface, 



In this section, we analyze the canonical problem of 
ICP in a leaky membrane illustrated in Fig. [3j A sym- 
metric binary electrolyte (D± = D, z± = ±z) passes 
from a reservoir of fixed concentration (c± = Co, <j) = 0) 
at x — through a weakly cation-selective leaky mem- 
brane with p s < through an ideal anion-blocking sur- 
face (F- = 0, zeA p F + = I, cf> = —V) at x = L, which 
could represent a non-leaky cation-exchange membrane 
or an electrode consuming cations, e.g. by electrodeposi- 
tion. We define I = A p J, as the total current that passes 
through the cross-sectional pore area A p , and solve for 
the transient current-voltage characteristics of the leaky 
membrane itself, not including interfacial polarization at 
either end. 



dine., ~. dd> ,„ ~, 
as well as 4>(0,t) = and 4>(l,t) = —V, where 

v= zeV 



kuT 



(18) 



(19) 



is the dimensionless applied voltage across the leaky 
membrane (not including interfacial voltage drops and 
series resistances [24j). 



Steady State for a Dilute Electrolyte 



Dimensionless Equations 

Following our prior work [THl |55] , it is convenient to 
transform the LMM for a symmetric, binary electrolyte, 



For constant diffusivity in a dilute electrolyte, the 
steady state can be solved analytically [26 . The con- 
centration profile is given by an implicit formula, 



C — ps ln(c) = 1 — Ix 



(20) 



D(c) 



z:z 

electrolyte 

reservoir 



+ + 
+ + 
+ + 

+ + 
+ + 
+ +1 
+ + 



= 



LEAKY MEMBRANE 



current / 



+ + - + 



+ 
+ 



c ± (x,t), 0(x,f) 
ze(c_-c + )= p s <0 







cation 
selective 
interface 

±+±±t± 
J+ttJ+ 

V + -V + - 

,;- + -v-t + - 

F_ = 

F + =I/zeA p 

= -V 



FIG. 3: Canonical problem for the Leaky Membrane Model 
in one dimension. Symmetric binary electrolyte transport 
from a reservoir through a weakly cation-selective leaky mem- 
brane to an ideally cation-selective surface, such as a (non- 
leaky) cation-exchange membrane or an electrode undergoing 
cation electrodeposition. Solutions appear in Figs. |4]|9] for 
steady and transient applied currents. 



with D — Tl 

boundary conditions yields 



Integrating by parts and applying the 



I = l-e- v D(e- v )-p s VD(e 



-V\ 



(c—ps lnc)— — dc. 
dc 

(23) 
If D varies significantly along the channel, the current- 
voltage relationship will begin to deviate from the ideal 
case. In particular, overlimiting current will increase if 
the diffusivity strongly increases with decreasing concen- 
tration, which is often the case. 

To illustrate effects of non-ideal thermodynamics via 
Eq. (|8| , we consider the Debye-Huckel theory of electro- 
static correlations in a dilute z : z electrolyte. The molar 
activity coefficient 7 of each ionic species is given by 



In 7 = 



M 2 

%TTskT\ D ' 



where e is the permittivity of the solvent and 



An = 



/ ekT 
2(ze) 2 c 



(24) 



(25) 



and current-voltage relationship 



1 



PsV 



(21) 



has the form of an equivalent circuit consisting of a 
diode, representing neutral-electrolyte concentration po- 
larization, in parallel with a shunt resistor, representing 
surface conduction. At low voltages and high bulk con- 
ductivity, the model describes the familiar linear Ohmic 
regime of bulk electrolyte transport, I ~ V(l — p s ). At 
high voltages and low bulk conductivity, however, the 
model predicts OLC, I ~ 1 — p s V, with a constant 
over-limiting conductance sustained by SC. This formula 
provides good fit of experimental data for conduction 
through a silica glass frit (a leaky membrane) from a 
reservoir to a Nation membrane in copper sulfate solu- 
tion [24] , although the over-limiting conductance also in- 
cludes strong effects of electro-osmotic flow [35] . 



Concentration-Dependent Diffusivity 

In most of our calculations, the diffusivity is treated 
as a constant, independent of concentration, which is 
strictly valid only for dilute solutions. In leaky mem- 
branes, however, significant concentration variations are 
possible, which alter the theoretical predictions. It turns 
out that the steady-state problem can still be solved 
exactly in implicit form. Re- writing Eq. (201 with a 
concentration-dependent diffusivity gives 



is the Debye screening length, assumed to be larger than 
the effective hydrated ion size. The activity is reduced by 
attractive interactions between an ion and its screening 
cloud as the ionic strength is increased, In 7 oc —^/c. The 
tracer diffusivity D® is taken to be constant, consistent 
with the moderately dilute range of validity for Debye- 
Huckel theory. 

d23l and (24), the effect of a non- ideal dif- 



D^_ £) ~ pa dM^ 

dx dx 



(22) 



Using Eqs. 

fusivity ([8]) can be found for varying initial ion concen- 
tration. In Figure [4] we see that for very dilute solutions 
(ImM) there is very little change in the current- voltage 
relationship and concentration profile. However, at larger 
concentrations (1M), there is a significant deviation. The 
Debye-Hiickel theory also breaks down and underpredicts 
the activity, so this example suffices to bound the trends. 
The overlimiting current is larger than expected from the 
ideal case and the depletion region is wider. This arises 
as a result of the diffusivity increasing in the depleted 
region leading to an increase in mass transfer. This cal- 
culation shows that it is generally necessary to account 
for thermodynamic corrections in the LMM for concen- 
trated electrolytes, although the qualitative results are 
similar with ideal solution theory, consistent with exper- 
imental data |2"4"1. 



Mathematics of Deionization Shocks 

Conductivity Dynamics at Constant Current 

At constant (over-limiting) current, a leaky membrane 
(or microchannel [20] ) supports the propagation of deion- 
ization shocks [!!]. Below, we will discuss how shocks 
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FIG. 4: Effects of non-ideal thermodynamics in the model 
problem of Fig. [3] Steady-state current-voltage relations 
(Top) and concentration profiles (Bottom) are shown for two 
initial ion concentrations, ImM (left) and 1M (right), using 
the concentration-dependent diffusivity D(c) from the Debye- 
Hiickel theory of electrolyte activity. 



arise in transient response of a finite system, but we first 
discuss the mathematics of shock propagation. It is con- 
venient to recast the equations in terms of the total ion 
concentration, k = c + + c_ , which is proportional to the 
total conductivity (for equal ionic diffusivities) , including 
both bulk and surface contributions. Note that this def- 
inition of k is different from the "bulk" conductivity Kb 
defined by Mani and Bazant [T5], set by the depletable 
co-ion concentration c = C-. 

In terms of the dimensionless conductivity k 
the LMM takes the form 



2c ' 



dk 



d 2 k 
dx 2 



p s I dk 

k? dx 



with boundary condition 



Ok 

dx 



(M) 



PsI 
«(M) 



= -I, 



(26) 



(27) 



The conductivity can be solved independently and the 
potential then obtained by integrating the current 



I =-k 



d'i 



(28) 



By rescaling Eq. (26) we can obtain a simple, quasilin- 
ear PDE 



dk d(k 1 ) d 2 k 
di dx dx 2 ' 



(29) 



for the rescaled conductivity k = k/y—p s I (assuming 



that p s is negative). This PDE resembles Burgers' equa- 
tion [HI [34], where k 2 /2 has been replaced by k _1 . While 
Burgers' equation can be transformed into the linear dif- 
fusion equation by the Cole-Hopf transformation [T51I16J . 
there does not seem to be a suitable linearizing transfor- 
mation for Eq. (29). However, the scaling used to achieve 



this simplification suggests that the surface charge den- 
sity may be as important as the applied current in affect- 
ing the resulting conductivity profile. In fact, it will be 
shown subsequently that p s significantly impacts the con- 
ductivity and voltage profiles, bringing about non-linear 
behavior. 

As with Burgers' equation, the long-time dynamics in 
an semi-infinite medium is dominated by nonlinear ad- 
vection (the second term in Eq. (29)), which leads 



to shock waves for most boundary and initial condi- 
tions. Diffusion plays only a secondary role in determin- 
ing shock structure [19] . Neglecting the diffusion term, 
we obtain a first-order quasilinear PDE, 



,dk 



dk 
d.r 







(30) 



of the same form as Amundson's model of chromatog- 
raphy, Eq. pi), which can be solved by the Method of 
Characteristics I31ITT1. 



For current in the +£ direction, nonlinear kinematic 
waves in the conductivity profile propagate in the —x 
direction. For characteristics originating at the end of 
the leaky membrane where conductivity variations are 
specified, k(0,t) = ko{t), the solution for x < is given 
by 



k = ko (i + k x) 



(31) 



For characteristics originating in the bulk material, the 
initial conductivity profile, k(x,0) = fci(x), evolves as 



«i \ x + 



t 



(32) 



These solutions are valid until characteristics nearly 
cross, as conductivity gradients become large and lead to 
shock waves. Multi- valued solutions ("wave breaking") 
for the conductivity profile are prevented by diffusion, 
which stabilizes the shock structure. 



Conductivity Shock Waves 

A dei onization shock is a traveling-wave solution of 
Eq. (29), k(x,i) = /(£) with £ = x — v s i, where v s is the 
shock velocity. Let k — /_oo and k = foo < /_oo be the 
conductivity asymptotes ahead (x — > — oo) and behind 
(x — > oo) the shock, respectively. The shock profile then 



satisfies the ordinary differential equation, 

- V.f + (r 1 )' = /", /(±00) = /±oo 



(33) 



Integrating once we obtain, 

-« s (/-/oc) + (/- 1 -/- 1 ) = / / (34) 

where we impose the boundary condition at £ = oo, 
where /' — > 0. If we also impose the boundary condition 
at £ = — oo, we obtain the shock velocity (a nonlinear 
eignevalue): 



pioneered by Amundson for transport problems): 



J oo J — 

Joo /— 



<0 



(35) 



for propagation directed toward the high conductivity re- 
gion, leaving behind a depleted region behind the shock. 
Equation (35 1 is the Rankine-Hugoniot jump condition 



expressing integrated mass conservation across the shock. 

As in the c(x,t) formulation [T!5], it is possible to in- 
tegrate (34) analytically to obtain the shock profile for 



k(x,i) in implicit form, but it is fairly complicated. A 
much simpler, approximate solution can be obtained by 
neglecting the nonlinear advection term ahead of the 
shock, 



m) 



\J— oo Joo)& 



for £ < 
for £ > 



(36) 



which corresponds to a truncated exponential profile 
moving at constant velocity, clearly seen in some of our 
simulations below (Fig. rol 1 = 5). This "diffusive wave" 
solution arises whenever an absorbing boundary propa- 
gates into a diffusing medium, as in dendritic electrode- 
position [57] , where metal deposition plays the role of sur- 
face conduction in rapidly removing cations from the bulk 
solution ahead of the wave. In the absence of flow, deion- 
ization shock waves are nonlinearly stable to conductivity 
perturbations [19] . due to a mathematical analogy with 
interface motion in diffusion-limited dissolution 1881. 



Transient Response to a Current Step 

A canonical problem of leaky membrane dynamics is 
the response to a current step for a dilute, symmetric 
binary electrolyte in a charge porous material. Three 



different dynamical regimes in the solution of Eqs. ( 15 ) 
( 19 1 can be identified: 



Zero Surface Charge: Neutral Electrolytes 

The limit of zero surface charge, p s = 0, corresponds 
to an "ideally leaky membrane" consisting of quasi- 
neutral electrolyte confined within the pores. The classi- 
cal Nernst-Planck equations apply, only with diffusivities 
corrected for the tortuosity and porosity. An exact series 
solution can be obtained by finite Fourier transform (as 



a = i - ix + 2 



oo 

£ 

n=0 






-Kt 



sinA„a:, (37) 



where X n = (n + |)tt. According to Sand [89"] , this clas- 
sical solution of the diffusion equation was first derived 
by H. F. Weber in 1879 [3D], who applied it to infer the 
diffusivity of Z11SO4 from the voltage transient after the 
interruption of steady current. The series can be trun- 
cated at a small number of terms without losing accu- 
racy for late times, i > 4- The series is non- uniformly 
convergent, however, and requires a diverging number of 
terms at early times. 

For early times or large currents, a more accurate and 
insightful similarity solution (which effectively sums the 
series) can be obtained by solving for jj| in a semi-infinite 
domain. After integrating ||| and applying the boundary 
conditions, the concentration profile is found to be: 



c = 1 + 21 v t r\ erfc ij 



1 



2 Vt 



(38) 



This famous result was first obtained by Sand in 
1901 [89] , who applied it to infer the diffusivity of CUSO4 
from observations of "Sand's time" [? ], isand, the time 
when the voltage diverges during constant over-limiting 
current. Solving Eq. (38) for c(l,i) = shows that in 
this case, 



^Sand — 



TV 

4J2' 



(39) 



At Sand's time, the concentration goes to zero at the 
selective surface, and the potential at that point is un- 
defined and corresponds to an infinite voltage. A crucial 
observation, however, is that this is voltage spike, a sig- 
nature of diffusion limitation, can be removed by surface 
conduction in a leaky membrane. 



Large Surface Charge: Ion Exchange Membranes 

The other extreme is the case of p s ^$> 1, which corre- 
sponds to a highly charged ion-exchange membrane with 
high electrochemical permselectivity. As shown in Fig- 
ure [5j high values of surface charge suppress large con- 
centration gradients, even under OLC. This lack of con- 
centration gradient results in a nearly constant potential 
across the system, meaning that under p s 3> 1 conditions 
the system is almost a purely controlled by the diffusion 
of the counter ions. Interestingly, the transient behav- 
ior of the case of very large p s behaves similarly to that 
of the case of p s = (as shown in the FFT solution), 
where transport is dominated by diffusion in the absence 
of a significant potential gradient. In a neutral medium 
(p s = 0) the ambipolar diffusivity (based on the diffusiv- 
ities of the counter and co-ions) determines the transient 
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FIG. 5: Salt concentration evolution after a current step in 
slightly leaky ion-exchange membranes (Fig. 3) with moder- 
ate (p a — — 1) and high (p s = —10) negative surface charge, 
just below (/ = 0.9) and far above (/ = 5) the limiting cur- 
rent. 



behavior, but the counter-ion diffusivity alone dominates 
at high charge density (p s » 1). 



Small Surface Charge: Leaky Membranes 

In a leaky membrane, the surface charge is relatively 
small, p s = 0(\) but plays an important role. Below 
the limiting current (I <C 1), a small dimensionless sur- 
face charge (0 < \p s \ <C 1) acts as a regular perturbation 
of the system, and the solution to our model problem 
remains close to the diffusive relaxation of a neutral elec- 
trolyte (37) for all times. Above the limiting current 



(I > 1), however, even a very small, but non-zero, sur- 
face charge (0 < \p B \ -C 1) acts as a singular perturbation 
that significantly alters the dynamics. The transient con- 
centration profile in our model problem for three different 
currents is shown in Figure [6]and several OLC voltage 
responses are given in Figure [7j Under OLC conditions, 
the ion concentration profile undergoes three stages: 1) 
Depletion, 2) Shock Propagation and 3) Relaxation. 

1. Salt Depletion. During this early stage the ion 
concentration is depleted at the selective surface 
and behaves similarly to the classic model (p s = 0) . 
This is more clearly shown in second half of Fig- 
ure [7] where time has been rescaled with respect 
to the classically derived Sand time, t Sand . With 
this rescaling it is clear that the large voltage in- 
crease, corresponding to full depletion, occurs at 
t\ = isand- The fact that time scales for the clas- 
sical case still apply when SC is taken into con- 
sideration is further shown by demonstrating the 




FIG. 6: Transient response of the salt concentration to a cur- 
rent step across a leaky membrane (p s — —0.01, Fig. 3). Just 
below the limiting current (a), the dynamics is dominated 
by linear diffusive relaxation, while nonlinearity begins to al- 
ter the concentration profile just above the limiting current 
(b). Well above the limiting current (c), the initial diffusion 
layer drives total salt depletion (stage 1), followed by a new 
dynamical regime of shock propagation (stage 2), ending in 
relaxation to steady state (stage 3). 



11 



1000 

900 
800 
700 
600 
'> 500 
400 
300 
200 
100 






/ 


/ / 




/ / 3 


- 


/ ' ' 
/ ' ' 
1 l / 

1 ' •' 


- 


/ / 2 / 


1 = 4 

1 = 6 

1 = 8 

1 = 10 




FIG. 7: Transient voltage of across a leaky membrane in re- 
sponse to a current step with increasing current (p s = —0.01). 
Top: Dimensionless voltage versus dimensionless time. Bot- 
tom: After rescaling time to Sand's time, three distinct dy- 
namical regimes (Fig. p| appear at high current: 1. depletion, 
2. shock propagation, and 3. relaxation to steady state. 



impact of p s in Figure [8] In this figure the voltage 
response is shown for decreasing surface charge. As 
the absolute value of p s decreases an order of mag- 
nitude the voltage drop increases about an order 
of magnitude. As the dimensionless surface charge 
continues to decrease the system grows closer to the 
classical result, as expected. 

2. Shock Propagation. After the co-ion concentra- 
tion is fully depleted at x = 1, a deionization shock 
appears and propagates away from the selective 
surface [20] with a stable exponential profile [19], 
given by Eq. (36). In the case of total salt deple- 
tion behind the shock (kqq = — p s -C 1) and unper- 
turbed conductivity ahead (k_oo = 1 — Ps > 1), the 
shock velocity, Eq. ( 35 1 , takes the simple form 



v. = -z =- (40) 

I- Ps 

which is equal to the (dimensionless) electromigra- 
tion velocity of co-ions [TH1 [57] , as required by mass 
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FIG. 8: The effect of varying the surface charge on the tran- 
sient voltage response to a current step (Fig. YA for I — 4. 



conservation across the shock, 
left behind. 



since no co-ions are 



Over the duration of stage 2 r 2 can be estimated 
as the time for the shock to move at velocity v s 
from the selective surface at x = 1 to the edge of 
the steady-state depletion region at x = I~ l . This 
implies the scaling: 



r 2 =(/- 1 -/- 2 



(1-P-) 



(41) 



Next we analyze the transient voltage during stage 
2. The dimensionless conductivity in the depleted 
region is approximately p s < 1, which dominates 
the total electrical resistance of the leaky mem- 
brane. The length of the depleted region at any 
time past ts&nd is equal to the shock velocity (/) 
times time. Therefore the resistance of the depleted 
region is I(t — ts a nd)/Ps- Using Ohm's law, the 
voltage thus scales as 



V 



\t — tSand) 
Ps 



(42) 



This scaling is verified at high currents in Figure [9j 
where Vp s /I is plotted against I(i—isand), leading 
to a data collapse of both stages 2 and 3 of the dy- 
namics. As the applied current increases, thereby 
strengthening the shock, the system closely obeys 
these scaling laws. At larger currents the depleted 
region nearly encompasses the entire system length 
with the shock propagating for a interval scaling as 
f 2 ~ I -1 after Sand time has been achieved. 

Aside: Shock Propagation at Constant Voltage. 
While the depiction region grows linearly with 
time under constant current conditions, it has been 
shown that depletion propagates as t 1 / 2 under 
constant- voltage conditions |2"3"] . This scaling can 
be easily revealed by estimating the shock as a mov- 
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applied currents. (p s = —0.01) 



FIG. 10: Sketch of the model problem flow aligned with the 
current, normal to the cation selective surface. Current is 
passed from a reservoir through a leaky membrane to a flow- 
through porous electrode that consumes cations, but allows 
neutral salt to pass by convection. 



ing step function, with a dep letion region length 
equal to 6(t). From Eq. (27 1 and the boundary 



conditions the voltage drop across the system is 
given as V = I(t) L \dx. Integrating over the step 
function gives V sw j[(l — <5)/k_oo + 6/iioo]- Since 
Koo "C ii-oc, V s=s —18/ps- The size of the depletion 
region is equal to the integral of the shock veloc- 
ity, or I = d5/di. Combining this with the voltage 
drop approximation, the depletion region length is 
found to propagate as 



5(t) ss y — 2p s Vt (constant voltage) (43) 

demonstrating the i 1 / 2 behavior. 



steady state by adding the times for all three stages 



f = (1 - MI' 1 + 



(45) 



In the limit of large currents, / 3> 1, the response 
time is dominated by the time for the shock to 
cross the full thickness of the leaky membrane (first 
term). 



Steady State with Normal Flow 



Relaxation to Steady State. Under moderate 
surface charge, SC only plays a dominant role when 
ion concentrations are very low, such as in the de- 
pletion region. Outside of this region, transport 
is dominated by linear diffusion. Once the shock is 
close to its final position (determined by the applied 
current or voltage) the concentration profile relaxes 
to the steady-state profile. As diffusion is the dom- 
inant transport mechanism, the transient scaling 
during this stage will be similar to the p s = case, 
solved earlier by FFT. The main difference is that 
the relevant length scale is not the total leaky mem- 
brane thickness, but rather the width of the steady- 
state diffusion layer I . As such the eigenvalues 
are rescaled to A n w (n 
scale for relaxation is 



|) 7T J, and thus the time 



T3 



1 



7T 2 / 2 



(44) 



In the previous section we examined how OLC cre- 
ates ICP by forcing ion depletion regions to develop. 
In this section we explore the effects of a uniform nor- 
mal flow u = U directed toward the selective surface on 
the depletion region and compute steady-state concen- 
tration profiles and current- volt age characteristics. This 
idealized situation shown in Fig |10| is relevant for flow- 
through porous electrodes [SS], as well as dominant nor- 
mal flow that exits through a small side outlet near the 
membrane [24] . Instead of a solid right wall, a perfect 
porous electrode is in place at x = L, which allows fluid 
and neutral salt to pass through while removing all excess 
cations. 



Concentration Profiles and Polarization Curves 



Finally, we are able to predict the total time to 



Starting from Eqs. d9])-(ll|, a pair of dimensionless 
equations for the steady state is achieved by averaging 
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over the cross section, 






Pe d£ 
dx 


d 2 c . d 2 $ 
dx 2 dx 2 


= 


d 
dx 


dx 



(46) 

(47) 



where the Peclet number Pe = UL/D controls the im- 
portance of convection relative to diffusion. Note that 
convection drops out of the charge balance (second equa- 
tion) because we assume a constant surface charge den- 
sity. For a flow-though electrode, the boundary condi- 
tions are more subtle. Only neutral salt can pass through 
the electrode, and therefore, to the right of x — L, uc + 
must equal wc_ . However, charge conservation within the 
membrane forces c+ > c_ between x — and x — L. As 
a result, a streaming current develops, I stream = —Up s , 
that accounts for this imbalance. The total current is the 
sum of the electro-diffusive fluxes, discussed above, and 
the streaming current: 



(C-Ps) 



dx 



Pep s 



(48) 



The anion electro-diffusive flux also vanishes, Eq. (181 



The equations can also be rearranged to determine the 
dimensionless total ion concentration (or conductivity), 
k — c — p s , instead of c: 



(I 



Pe— 

dx 

Peps) 
«(0) 



d?k 

dx 2 

dk 

dx 

1-/5 



(1) 



p s (I + Pep s ) dk 
k 2 dx 

p s (I + Pep s ) 



«(1) 



(49) 

(50) 
(51) 



Once this boundary value problem is solved for k(x), the 
potential profile and voltage are obtained from the cur- 
rent relation by a simple integration: 



(x) = -(I + Pep s 



ds 

k(s) 



(52) 



With flow, the analytical solution becomes more chal- 
lenging, so concentration profiles and current-voltage re- 
lationships are found numerically. 



In Figure 11 several concentration profiles are shown 
for varying voltage and Pe values. As shown 

previously [26], increasing the applied voltage increases 
the amount of depletion. The addition of convection 
pushes all ions toward the membrane, and the linear 
concentration profile of the quasi-steady diffusion layer 
gives way to the exponential profile of a diffusive wave, 
propagating against the flow [115] [57] . As the flow rate 
increases the concentration profile appears more shock- 
like with a decreasing shock width. Additionally, as the 
flow rate increases, the depletion region shrinks, requir- 
ing a higher applied voltage to maintain its size. At high 
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FIG. 11: Concentration profiles in a leaky membrane (p s = 
—0.01) with normal flow (Fig. 10 1 while varying (a) the volt- 



age with a high flow rate, Pe = 5, or (b) the flow rate at high 
voltage, V = 40. 



Pe the steady-state concentration profile converges to the 
propagating shock solution, Eq. (36), where the uniform 
flow holds the shock in place. 



Figure [12] shows the current-voltage relationship for a 
case of strong flow, Pe = 5. The shape of the curve is 
noticeable different from the case with no convection [261. 



Eq. (21), with a delayed, curved transition to the over- 



limiting regime. However, at higher voltages the over- 
limiting current eventually becomes linear, similar to the 
no-flow case. 



Energy Cost of Deionization 

In the regime of over-limiting current, the flow-through 
electrode is continuously dcionizing the fluid as it passes 
through the leaky membrane. This setup could have ap- 
plications to flow-through capacitive desalination [55], 
where the flow channel is filled with a porous medium 
or microchannel array, and it is a simple first approxi- 
mation for the cross-flow geometry of shock electrodial- 
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over-limiting current as the surface charge is increased. The 
vertical shift of the curves is due to streaming current from 
the convection of counter-ions screening the pore charge of 
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ysis [M] discussed below. For these applications, the 
model provides a simple case to analyze the energy cost 
of de-ionization. 

The energy per volume of initial solution processed, 
E v , is equal to input electrical power (IV) divided by 
the volumetric flow rate (Q): 



IV IV 

Ev - — -2kTc Q — 



which has the dimensionless form, 



E,, 



2kTc 



IV 
Pe" 



(53) 



(54) 



The energy cost, E v , is a function of the surface charge 
density, p s , the applied current, /, and the velocity, Pe, 
each in a suitable dimensionless form. 

Model predict ions for a fixed system geometry are 
shown in Figure 
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In these plots E v is shown versus 
varying Pe and I for two different values of p s (-0.01 and 
-0.0001). A black line indicating when the depletion re- 
gion has formed (c = 0.001) is placed on top of these 
surface plots. Below this line c is less than 0.001 and the 
energy per volume increases. The energy profile appears 
very similar in these two plots, with the values differing 
by a factor of 100. In the ID case, surface conduction 
goes as p s V. Therefore as p s decreases by a factor of 100, 
the voltage must increase by a factor of 100 to maintain 
the same level of conduction. This increase in necessary 
voltage is what leads to the 100-fold increase in energy. 
In order to reduce the energy per volume required for 
deionization, the depletion region should be as small as 
possible. This is the reason that increasing the flow rate 
(Pe) decreases E v . Similarly, increasing the applied cur- 




FIG. 13: Energy cost of deionization for the ID model of 
normal flow through a leaky membrane and porous electrode. 
The black line indicates where a depletion region has formed. 
Below this line the outlet concentration, c, is 0.001 or less, 
(a) p s = -0.01; (b) p s = -0.0001 



rent past the point of early depletion wastes energy and 
increases E v . In this ID model with uniform flow, the 
fluid recovery fraction, or ratio of deionized to incoming 
fluid volumes, is 100%. High fluid recovery is a hall- 
mark of flow-through separations in porous media, but 
the shock phenomenon provides an opportunity for effi- 
cient separations in cross flow, perpendicular to the cur- 
rent, which we analyze next with a 2D model. 



Steady State with Cross Flow 

Fractionation by Deionization Shocks 

The formation of a deionization shock represents a dy- 
namic, "membraneless" separation of salty and deionized 
solution, which can be exploited for water purification, 
brine concentration, or other separations by fractiona- 
tion in cross flow. In contrast to traditional electrodialy- 
sis (ED), in "shock electrodialysis" [23] there is no fixed 
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FIG. 14: Sketch of one layer of a cross-flow shock electrodial- 
ysis system for water deionization and brine concentration. 
Current in the normal (vertical) direction passes though a 
negatively charged porous material (leaky membrane) sand- 
wiched between two cation-exchange membranes (or other 
cation selective layers). Water flows in the perpendicular di- 
rection and is split into brine and deionized streams upon 
exiting. 



physical, chemical or electrostatic barrier between the 
two regions that spontaneously form in a homogeneous 
porous medium. The strong localization of the salt con- 
centration jump in the shock and its ability to propagate 
to a desired position enables separation in cross flow. 



The basic idea, sketched in Figure 14 is to drive fluid 
flow through the leaky membrane in one direction, while 
driving over-limiting current in the perpendicular direc- 
tion between ion perm-selective membranes. The deion- 
ization shock propagates in the cross flow to form a 
boundary layer of strong depletion, which extends across 
a fresh water collection outlet on the downstream end 
of the leaky membrane. If the leaky membrane is sand- 
wiched between identical ion-exchange membranes, then 
an enrichment diffusion layer also forms on the other side, 
which is collected in a brine stream, separated from the 
fresh stream by splitting the flow leaving the leaky mem- 
brane. This layered structure is the basic building block 
of a scalable shock ED system with electrode streams on 
the ends to sustain the current. 

Several parameters will affect the efficacy and the ef- 
ficiency of such a device, including the surface charge of 
the leaky membrane, the geometry of the device, and 
the applied current or potential. In order to understand 
how these parameters relate to each other, we use a 2D 
LMM. Analytical solutions with nested boundary layer 
structure are possible in the relevant regime of fast cross 
flow, where the deionization shock and enrichment re- 
gions remain well separated |91j . but here we focus on 
numerical solutions for finite geometries in a wide range 
of operating conditions. The purpose of this model is to 
understand in a simple way how forced convection and 
SC affect ion transport while providing design guidelines 
to optimize the system. 



Two Dimensional Model 

Consider a leaky membrane of height, h, and length, 
L, where x S [0,h] and y € [0, L\. Uniform flow with 
velocity, u, in the y direction originates from y — 0. In 
this 2D model, we neglect axial diffusion, which is valid 
beyond a distance D/u from the inlet, which is small if 
Pe y = uL/D ^> 1. In this regime, convection dominates 
in the y-direction, and diffusion in the x-direction. The 
LMM conservation equations then take the form 



dc + 
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oy 



D 



= D 



d 2 c+ 

dx 2 
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Nondimcnsionalizing as before, with x = f 
the dimensionless conservation equations are 



(55) 



and y = I, 



uh 2 dc d 2 c 
LD dy dx 2 



d_ 

dx 



(5- Pa) 



d<p 
dx 



Ps dx 2 ' 

= 0. 



(56) 

(57) 



In our previous examples (Figs. [3J [10| , the leaky 
membrane was in contact with a reservoir of constant 
concentration at one end (x = 0) and a cation-selective 
membrane or porous electrode at the other end (x = h). 
In this case, anions are blocked at both ends (x = and 
X = h), which implies Neumann type conditions, 



x = 

x = l 
y = 



-V 



<91nc 



= 0, 
5=1. 



dx 
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dx 



dx' 
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(58) 

(59) 
(60) 



with Dirichlet boundary conditions for the concentration 
at the inlet and potential at the membranes. The current 
density (per area) is no longer uniform, 



J(y) 



dx' 



(61) 



where J = „_{k 



'z Dc 1S ^ ne dimensionless current density in 
the x-direction. J must be integrated over the membrane 
area to obtain the total current 



J(y)dy 



(62) 



Noticing that several parameters of interest are lumped 
together, the conservation equation can be rewritten in 
terms of the Peclet number 



Pe 



uh 

~D 



(63) 
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and a new axial length variable, 



II 



yD_ 

uh 2 



hPe 



(64) 



scaled to the entrance length for the convection-diffusion 
boundary layer, uti 2 / D, as usual in the analysis of forced 
convection in a channel or pipe (28] . With this change of 
variables, 



dc 

Of) 



dx 2 



d 2 (f> 
' dx 2 



(65) 



the conservation equation becomes the same as the ID, 
transient equation, Eq.( 15). Therefore, the solutions 



from the previous section can be reworked and applied 

here. 

Example concentration profiles are shown in Figure |l~5"| 
for p s — —0.01,-0.05 and V = 30. As expected, in- 
creasing the surface charge increases the size of the de- 
pleted region, S. Here 5 was taken to be the point where 
c = 0.001. It is also important to note that at around 
y = 0.1 the concentration has reached its steady state 
value and no further depletion occurs. By setting y = 1 
these plots can be used to examine the outlet concentra- 
tion distribution. For example, in the case of p s = —0.05 
(Figure 15 o) the outlet can be fractionated at x = 0.25 
and if y > 0.1 only depleted fluid will be collected. This 
analysis can be used to determine the best geometry and 
flow rate. In order to maximize the flow rate, y should be 
minimized. In order to get full depletion, y should exceed 
the dimensionless distance to achieve steady state, which 
for these two cases is around 0.1; in other words, full 
depletion occurs at roughly 10% of the entrance length. 
Additionally, scaling up the system will not be a linear 
process, since y ex A. 




FIG. 15: Steady concentration profile in a simple 2D model 
of the shock electrodialysis device of Fig. [14] at high voltage 
V = 30 for varying dimensionless surface charge in the leaky 
membrane: (a) p„ = —0.01, b) p s — —0.05. 



Design Principles for Shock Electrodialysis 

An important characteristic of this extraction tech- 
nique is the energy required to create the necessary de- 
pletion region. For a volumetric flow rate of Q, with cur- 
rent density, I, and voltage, V, the energy E v required 
per volume of depleted solution is given by 



E„ — 



IV _ 2k B Tc IV 



(66) 



The parameters I, V, and 5 are related through the con- 
servation and current equations. If y is far enough down- 
stream to reach steady state, then V and <5 can be found 
based solely on / and p s . As a result, it is useful to 
consider the dimensionless energy efficiency, 



E„ 



E v Pe IV 
2kgTcQ 8 



(67) 



A plot of E v over a range of V and p s values is shown in 
Figure |16| This plot was generated for a system with a 
maximum length of y = 0.1, corresponding to about 10% 
of the entrance length. At lower to moderate applied 
voltages, increases in the surface charge density lead to 
decreases in the depletion energy. While increases in \p s \ 
will lead to increases in / the corresponding increases in 
5 are sufficient to lower the required energy. However, 
at higher applied voltages the balance is shifted and the 
increase in power cost overwhelms the efficiency gained 
by creating a larger depletion region. Once V and p s have 
been determined, the energy efficiency can be calculated. 
This efficiency can be enhanced by properly designing 
the system geometry. For example, the larger the aspect 
ratio, L/h, the lower the energy efficiency, as seen in Eq. 
p). 



In addition to energy requirements, in order to develop 
a practical device, the volume of depleted fluid relative 
to the incoming fluid (the "water recovery" percentage 
in desalination) must be considered. A very efficient de- 
vice that only depletes 1% of an incoming stream may 
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FIG. 16: Dimensionless energy E v (color contours) per vol- 
ume of de-ionized fluid in the model of Fig. |15| versus the 
dimensionless surface charge p s and dimensionless voltage V, 
evaluated at y = 0.1. 
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FIG. 17: Water recovery S (ratio of deionized to incoming 
fluid volumes) in the model of Fig. 15 which increases with 
increasing V and p 3 . 



not be particularly desirable. In this case many passes 
would be required to achieve a sufficient amount of de- 
pleted solution. Recovery in this model corresponds to 
the size of the depletion region. 5. A plot of S versus 
V and p s is shown in Figure |l7[ under the same con- 

The region of highest recovery 



ditions as in Figure 16 



does not correspond to the region of lowest energy per 
volume (Figure 16). Therefore a balance must be struck 



between the two values, depending on the requirements 
of a desired system. 

In order to fully design a practical device that utilized 
SC in flow, one other parameter needs to be addressed. 
Throughout this study the parameter p s has played an 
important role. However, this parameter is a function of 
the initial anion concentration, Cq. The power of SC goes 
up as Co goes down. In order to have an effect on higher 
ion concentrations, the volume surface charge density, p s 



of the material should be increased. This can be done 
by either altering the surface charge of the material, a s 
or decreasing the pore size. For instance, a typical sil- 
ica bead in water has a surface charge density of about 
-0.001 coul/m 2 92j. In a ImM solution, a porous struc- 
ture of these beads with a pore size of 10 /im will result 
in p s — —0.001. However, if the pore size decreases to 
100 nm, then p s = —0.1 and SC plays a more dominant 
role. Based on this analysis, the smaller the pores, the 
better the de-ionization. However, the energy analysis 
conducted here did not take into consideration the force 
to pump the fluid through the porous material. As the 
pore size decreases the pump energy required increases. 
As a result, decreasing the pore size may not be the best 
solution. Alternatively, the surface of the porous ma- 
terial can be altered to create a more negative surface. 
In order to maximize the energy efficiency, the device 
should be designed with a high aspect ratio. Addition- 
ally the velocity should be maximized such that y de- 
fined in Eq. (64) is kept low but at a steady state value. 



Based on Figures [16] and [TTJ the applied current should 
be above the limiting value but low enough such that the 
energy per volume and water recovery are at acceptable 
levels. In this manner, an efficient SC-flow device can be 
developed using simple materials. 



Conclusion 

Unlike ideal ion-exchange membranes, which maintain 
a large conductivity of counter-ions, the conductivity of 
"leaky membranes" with larger pores and/or smaller sur- 
face charge densities can vary significantly in response to 
a large applied voltage. The surface conductivity, which 
remains even if the bulk salt is depleted, provides a mech- 
anism for over-limiting current, faster than diffusion. 
This can lead to a macroscopic region of salt depletion 
behind a propagating deionization shock, which opens 
new possibilities for nonlinear electrokinetic separations 
in porous media. Building on recent work [2U1 EH] > 
we formulate a general Leaky Membrane Model and de- 
rive representative analytical and numerical solutions for 
finite domains. We focus on the simplest situation of 
a symmetric binary electrolyte in a leaky membrane of 
constant surface charge density, uniform pressure-driven 
flow, negligible hydrodynamic dispersion, and no electro- 
osmotic flow. 

Relaxing these assumptions and deriving suitable mod- 
ifications of the model provide challenging avenues for 
research. For example, charge regulation in a multi- 
component electrolyte due to specific adsorption of ions 
is a classical source of nonlinearity [3l 161] . which in 
leaky membrane can lead to over-limiting current by 
"current-induced membrane discharge" [15] ■ The LMM 
with charge regulation could have relevance for electroki- 
netic remediation in soils [30H33] . as well as ion trans- 
port in biological cells. The dynamics of charged colloids 
in leaky membranes may also lead to interesting nonlin- 
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ear dynamics, generalizing shock waves in capillary elec- 
trophoresis [34-38J. The LMM may also improve the 
accuracy of porous electrode theories, which currently 
assume electroneutrality in the solution phase and ne- 
glect surface conduction [5T], which already account for 
capacitive charging of double layers [93] with Faradaic 



reactions [7H [75] and specific adsorption of intercalation 
reactions [65j [66] , but generally neglect surface conduc- 
tion. In all of these situations, perhaps the most difficult 
and important extension of the LMM will be to account 
for electro-osmotic flow and associated dispersion phe- 
nomena at the macroscopic scale [26] [53] . 
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